Demographic responses to climate‐driven variation in habitat quality across the annual cycle of a migratory bird species

Abstract The demography and dynamics of migratory bird populations depend on patterns of movement and habitat quality across the annual cycle. We leveraged archival GPS‐tagging data, climate data, remote‐sensed vegetation data, and bird‐banding data to better understand the dynamics of black‐headed grosbeak (Pheucticus melanocephalus) populations in two breeding regions, the coast and Central Valley of California (Coastal California) and the Sierra Nevada mountain range (Sierra Nevada), over 28 years (1992–2019). Drought conditions across the annual cycle and rainfall timing on the molting grounds influenced seasonal habitat characteristics, including vegetation greenness and phenology (maturity dates). We developed a novel integrated population model with population state informed by adult capture data, recruitment rates informed by age‐specific capture data and climate covariates, and survival rates informed by adult capture–mark–recapture data and climate covariates. Population size was relatively variable among years for Coastal California, where numbers of recruits and survivors were positively correlated, and years of population increase were largely driven by recruitment. In the Sierra Nevada, population size was more consistent and showed stronger evidence of population regulation (numbers of recruits and survivors negatively correlated). Neither region showed evidence of long‐term population trend. We found only weak support for most climate–demographic rate relationships. However, recruitment rates for the Coastal California region were higher when rainfall was relatively early on the molting grounds and when wintering grounds were relatively cool and wet. We suggest that our approach of integrating movement, climate, and demographic data within a novel modeling framework can provide a useful method for better understanding the dynamics of broadly distributed migratory species.


| INTRODUC TI ON
The dynamics and trends of migratory animal populations depend on environments encountered across the annual cycle (Hostetler et al., 2015). For this reason, identifying networks of migratory connections represents a critical first step to the conservation of these species (Ruegg et al., 2020). Rapid advances in tracking technology have allowed for identification of these connections for a growing list of migratory animal species (McKinnon et al., 2013;Ruegg et al., 2014;Rushing et al., 2014). However, the value of this information for conservation is relatively limited without also identifying key environmental stressors at different points of the life cycle to better understand the consequences of migratory connections for population dynamics (Mancuso et al., 2021;Saracco & Rubenstein, 2020;Wilson et al., 2018).
Studies aimed at parsing the relative importance of demographic and environmental drivers of population change across the annual cycle have been based on a model of two sedentary life history stages (breeding vs. non-breeding) separated by brief spring and fall migratory periods (Rushing et al., 2017;Saracco & Rubenstein, 2020;Woodworth et al., 2017). However, additional periods of movement during the annual cycle of migratory birds may be more common than previously recognized (Pyle et al., 2018;Ruiz-Gutierrez et al., 2016), and these complexities should also be considered in full annual cycle models (Pageau et al., 2020;Pyle et al., 2018). For example, black-headed grosbeak (Pheucticus melanocephalus; Figure 1) is 1 of at least 19 bird species that breeds in seasonally arid regions of the western United States and has been documented migrating to the North American monsoon region (southwestern United States and northwestern Mexico) in late summer to undertake its annual prebasic molt prior to proceeding to overwintering areas farther south (Pageau et al., 2020;Pyle et al., 2009;Rohwer et al., 2005;Siegel et al., 2016).
Black-headed grosbeak populations have increased across much of their breeding range over the past 50 years; however, populations in California and other parts of the southwestern United States have tended to decline (Sauer et al., 2020). We suggest that spatial and temporal variation in black-headed grosbeak population dynamics and trends may reflect climate-driven patterns of vegetation phenology and productivity at breeding, molting, and wintering sites.
Identifying effects of climate variation on black-headed grosbeak populations and their habitats will be critical for the conservation of this and other species that utilize seasonally arid habitats across western North America given observed and predicted climatic shifts due to climate change. Some observed and predicted effects of climate change in areas used by this species across the annual cycle include severe drought and reduced vegetation productivity on breeding grounds (Diffenbaugh et al., 2015;Goulden & Bales, 2019;Trujillo et al., 2012;Williams et al., 2020), changes in the timing and amount of rainfall on the molting grounds (Cook & Seager, 2013;Grantz et al., 2007;Méndez-Barroso et al., 2009;Pascale et al., 2017), and warmer drier conditions on the wintering grounds (Karmalkar et al., 2011;Neelin et al., 2006).  (Wang et al., 2016), with a goal of assessing potential environmental drivers of vital rates and population dynamics in two breeding regions, the Sierra Nevada mountain range of California (hereafter Sierra Nevada) and along the coast and Central Valley of (hereafter Coastal California). The IPM framework allowed us to estimate and model recruitment rates using age-specific capture data from adult birds, and adult survival rates based on adult capture-mark-recapture data for both regions. We assessed effects of the following climate covariates on grosbeak adult survival and recruitment: (1) drought on breeding grounds, (2, 3) drought and rainfall timing on molting grounds, and (4) drought on wintering grounds. We delineated molting and wintering grounds based on locations and habitat relationships of four GPS-tagged birds that were captured and subsequently recaptured at MAPS stations.

T A X O N O M Y C L A S S I F I C A T I O N
year. All stations were operated during the breeding season (May 1 to Aug 9). Captured birds were aged as adults (after hatching year) or juveniles (hatching year). To the extent possible, adults were further microaged as either yearlings (second year birds) or older adults (after second year birds) based on criteria in Pyle (1997; Figure S1).
Our analysis included 880 year-unique captures of adult birds in the Sierra Nevada region and 2387 year-unique captures of adult birds in the Coastal California region. Of adult birds, 70% were microaged as either yearlings or older adult birds.

| Climate and vegetation data
To characterize annual climate variation and to link climate to demographic rates, we extracted point-level data from the ClimateNA database ver. 7.00 (Wang et al., 2016) for the 44 breeding sites (i.e., MAPS stations) and for 1000 random points in each of the non-breeding ranges. To control for spatial variation in climate covariates related to geographic gradients and habitat, we used differences between annual values and mean values for each point for the 1971-2000 normal period (referred to as deviations, hereafter). We selected climate covariates based on our preconceived ideas that water availability would be a strong predictor of annual variation in habitat quality and its effects on avian demographics. To characterize drought conditions on breeding, molting, and wintering grounds, we used Hargreave's climatic moisture deficit (CMD; Hargreaves & Samani, 1982), calculated as the difference between atmospheric evaporative demand and precipitation. Note that more positive CMD values indicate more severe drought. For the breeding grounds, we used summed CMD values from the overwintering F I G U R E 2 Distribution of black-headed grosbeaks that breed in two U. S. Bird Conservation Regions (Coastal California [light green] and the Sierra Nevada [darker green]) across the annual cycle. GPS-tagged birds were recaptured at sites in Yosemite National Park (YOSE; three birds, two MAPS stations) and Golden Gate National Recreation Area (GGNRA; one bird). Locations of these birds during the non-breeding season were separated into two distinct periods. Molting locations of individuals (more than one if they moved within-season) are indicated by dark orange shapes (circles for Yosemite, squares for Marin) within an orange molting range defined as a shrub and forest habitat-filtered minimum convex polygon +0.5-degree buffer surrounding August through mid-October locations. Overwintering locations are indicated by blue shapes (per above, based on region, with only three birds with winter data) within a similarly defined blue region surrounding late October through March locations (Artwork by L. Helton.) period (Oct-Apr) prior to the grosbeaks' arrival on breeding grounds, which is the period when most annual precipitation occurs in these two California regions; for the molting grounds, we used summed CMD values from Apr to Sep, encompassing a 6-month time window spanning a period prior to and including typical monsoon months when black-headed grosbeaks were present in that region; for wintering grounds, we summed CMD from April of the previous year to March of the current year. For molting grounds, we also calculated a ratio of early season (Jun and Jul) to late season (Aug-Sep) rainfall, to broadly characterize monsoon phenology.
To understand linkages between climate and vegetation, we extracted vegetation phenology and greenness data from the Land Cover Dynamics MCD12Q2 v.6 data product Gray et al., 2019), which is based on 0.5-km-resolution Moderate Resolution Imaging Spectroradiometer (MODIS) data from the NASA Terra satellite (http://terra.nasa.gov/). The MCD12Q2 v.6 data product uses the two-band Enhanced Vegetation Index (EVI) as the basis of phenology and greenness metrics; EVI is a composite metric that incorporates structural and seasonal components of vegetation greenness (Glenn et al., 2008;Potithepa et al., 2010). We considered two aspects of vegetation dynamics to encapsulate overall greenness and phenology. To quantify greenness, we used the integrated 8-day interval EVI over the main growing season (EVI area). To measure phenology, we used EVI maturity date, defined as the estimated date at which 90% of the EVI amplitude was attained. We derived site-specific annual deviations from mean values for all covariate values to account for spatial variation in overall timing and greenness. We fit linear models between mean EVI area and maturity responses and each mean climate covariate value, assuming normal responses, to assess climate effects on vegetation greenness and phenology.

| Integrated population models
We developed a Bayesian integrated population model (IPM) based on models of age-specific capture data and capture-mark-recapture data from MAPS stations and climate covariates. As a first step, we estimated annual indices of adult abundance for the two breeding regions based on a model of adult captures, denoted as AHY j,k,t . We assumed adult captures at station j, region k, and year t were distributed as AHY j,k,t ∼ NegBin r, p j,k,t , where r and p j,k,t relate to the mean count, AHY[j,k,t] , based on p j,k,t = r∕r + AHY [j,k,t] . We then modeled the mean parameter, AHY[j,k,t] , with the generalized linear mixed model: where 0[k] are fixed region-specific intercepts; 1 is the regression coefficient for an effort ef j,t effect (summed net-hours relative to average net-hours for the station); yr k,t is a random region × year effect distributed as Norm 0, 2 k ; and sta j is a random station effect distributed as Norm 0, 2 . We then derived the adult abundance index as â d k,t = e 0[k] +yr k,t , which was used to inform the population state of the IPM based on means and standard deviations of posterior MCMC samples: where the n k,t represents the population state variable of the IPM, derived as the sum of the numbers of survivors, nsurv k,t , and the numbers of new recruits, nrecr k,t . For the initial year (t = 1992), we assigned vague prior distributions for nsurv k,1992 and nrecr k,1992 . (2020). We assumed the probability of adult bird i at station j and region k in year t being a yearling bird to be distributed as a Bernoulli random variable, Y i,j,k,t ∼ Bern pY i,j,k,t , with the probability parameter modeled with the logit-linear model: logit pY i,j,k,t = 0[k] + yr k,t + sta j .
We derived estimates of the proportions of yearling recruits for each region and year based on back-transforming the region × year components of the model: p Y k,t = e 0[k] +yr k,t ∕ 1 + e 0[k] +yr k,t and estimates of the numbers of new yearling recruits as â d1 k,t =pY k,t ×âd k,t and numbers of surviving adults as â d2 + k,t =âd k,t −âd1 k,t . These estimates informed the recruitment and survival components of the population dynamics model: and Additionally, we derived an observed recruitment variable, ĝ k,t =â d1 k,t+1 ∕âd k,t , which we also assumed its posterior mean to be to inform the recruitment rate parameter.
The survival probability parameter of the population dynamics model, k,t , was informed by a state-space version of the Cormack-Jolly-Seber (CJS) model that accounts for transients (i.e., individuals that permanently emigrate after initial capture; Saracco et al., 2012) applied to CMR data from i = 1, …, M adult birds captured at MAPS stations. We modeled covariate effects based on generalized linear mixed models:  (ELR_m). We used CMD_m_res, rather than the original CMD deviation variable for the molting grounds, as it was moderately correlated with CMD_w (r = .62) and an initial analysis of the survival model suggested a stronger relationship for winter drought than for molting grounds drought. Thus, CMD_m_res can be interpreted as drought effects on molting grounds independent of drought effects common to both winter and molting grounds (Dormann et al., 2013;Graham, 2003).
Additional parameters in the CJS models included residency probability, i,j,k,t , and parameters describing the observation process, capture probability, p i,j,k,t , and probability of predetermining a newly banded bird as a resident bird, i,j,k,t (i.e., recapturing a bird ≥7 days apart in the season it was banded; Saracco et al., 2012). For residency probability, we defined a logit-linear model with random region × year effects; for the observation parameters, we defined logit-linear models that included intercepts and zero-mean random station effects.

| Movement and connectivity
All four GPS-tagged birds for which we recovered data were males and showed similar between-season movements ( Figure 2) Figure S1). While the range of dates spanned and resolution of the archived location data was greater for the two yearling birds (5-10 days) compared to the older birds (16-40 days), both younger birds made large-scale movements across a time period during which the older birds appeared to be largely stationary (Sep 5 to Oct 15).

| Climate variation and relationship to vegetation dynamics
Oct-Apr climatic moisture deficit deviation on the breeding grounds was highly variable among years with drought conditions tending to

| Population dynamics and demography
Adult apparent survival probabilities varied substantially among years, particularly for the Sierra Nevada populations (Figure 5a).
Recruitment rates were also annually variable for Coastal California but relatively stable across years for the Sierra Nevada (Figure 5b).
Average demographic rates were similar between the two regions  (1992-1993 and 1994-1995).

| Climate-demographic rate relationships
We found only weak support for most demographic rate-climate relationships that we considered (Table 2). Effect sizes, in several cases, were of a magnitude that could have large consequences on population size; however, precision of estimates was low, precluding strong inferences on climate effects on survival. We found strongest evidence of covariate effects with respect to recruitment rates for the Coastal California populations. Recruitment in this breeding region was higher following relatively cool and wet years on the wintering grounds (lower CMD_w; Figure 6a) and years with relatively more rainfall early in the Monsoon season on the molting grounds (ELR_m; Figure 6b). Our study provides a model of how tracking, environmental, and demographic data can inform temporal variation in habitat quality and patterns of population change at regional scales. Although based on may more accurately be considered to be protracted migrations than within-season movements. Support for non-breeding season movements has been reported for other migratory landbird species (Cormier et al., 2013;Mancuso et al., 2021;Ruiz-Gutierrez et al., 2016;Stutchbury et al., 2016), and flexibility to move to new sites within seasons may be a critical adaptation for tracking uncertain environmental conditions. Both yearling males that showed withinmolting season movements initially settled in habitats that were less green than those settled in by the older adults. It is possible that older males exclude younger males from preferred greener habitat.
However, apparent age differences could reflect alternative habitat selection strategies with younger birds selecting multiple habitats that green up at different rates ( Figure S2). Patterns of habitat use may also vary among years with habitat breadth expanding or contracting depending on conditions (Pyle et al., 2009). Finally, apparent age differences could also simply reflect lower availability of greener habitats in the year the yearlings were tracked (2017) compared to the year the older adults were tracked (2015). Clearly, additional data across age and sex classes across multiple years will be needed to fully elucidate non-breeding habitat needs.
We found clear signals of climate variation on habitat, with greenness and maturity date dependent on drought conditions at sites across the annual cycle and on rainfall timing at sites on molting grounds. In all seasons, more severe drought was associated with less greenness. Trends in drought across the black-headed grosbeak breeding range suggest that more severe and multi-year drought years are becoming more common, which will likely have implications for habitat quality (Breshears et al., 2005;Diffenbaugh et al., 2015;Goulden & Bales, 2019;Trujillo et al., 2012;Williams et al., 2020). Monsoon precipitation on molting grounds may also be declining (Pascale et al., 2017) and becoming more spatially (Demaria et al., 2019) and temporally (Meyer & Jin, 2017) variable, which could then affect spatiotemporal patterns in habitat variation (Méndez-Barroso et al., 2009).
Drought conditions were also associated with vegetation phenology. More severe drought on the breeding grounds was associated with earlier maturing vegetation, while more severe drought on molting and wintering grounds was associated with later maturing vegetation. A trend toward earlier maturing vegetation on breeding grounds could yield mismatches in the migration and breeding phenology of birds (Mayor et al., 2017); however, populations may also adapt to climate trends in space and/or time, and in some circumstances demographic rates can improve under warmer, drier conditions (Saracco et al., 2019;Socolar et al., 2017). We also found evidence that rainfall timing on molting grounds influences greenness and vegetation phenology. Although there was no trend evident in the our rainfall timing metric over the years we considered here, other studies have suggested a trend toward later monsoon rainfall (Cook & Seager, 2013;Grantz et al., 2007)  times of year where water is a limiting resource can adversely affect vital rates of migratory birds (Dugger et al., 2004;Rockwell et al., 2017;Saracco et al., 2018;Sillett et al., 2000). Thus, we expected drought effects on demographics might be greatest on the wintering grounds as the dry season peaks late in winter prior to spring migration and survival of young and adults may consequently suffer. While point estimates of effects of wintering grounds drought on adult survival were in the expected direction (i.e., negative effects of drought) for both breeding regions, credible intervals were broad, limiting inferences. Drought on the wintering grounds was more strongly correlated with declines in recruitment, although this relationship was limited to the Coastal California populations, which showed greater variation in recruitment and years of large population growth corresponding to high recruitment rates.
Recruitment in the Coastal California region was also associated with rainfall timing, with years having relatively more rainfall early in the Monsoon season tending to have higher recruitment than in years tending to have later rainfall. Given that earlier rainfall was associated with later vegetation maturation on molting grounds, it could be that such conditions were differentially favorable for first-year survival in such years, as young birds tend to migrate and molt later than adults (Ortega & Hill, 2020). Weak and/or sometimes conflicting relationships between demographic parameters and breeding and molting grounds climate covariates may reflect water not be a limiting resource at these life cycle phases. For example, productivity may decline for some bird species in coastal California in drought years (Chase et al., 2005); however, blackheaded grosbeaks finish breeding relatively early in that region F I G U R E 6 Relationships between recruitment rate in the Coastal California region and drought conditions on the wintering range (a) and ratio of early-to-late season rainfall on the molting range (b) and so may not be negatively impacted by drought, which may have stronger effects on habitats later in the season. Indeed, nest survival was higher for black-headed grosbeaks in relatively warm years at a site on the edge of the eastern Sierra Nevada/western Great Basin (Becker & Weisberg, 2015), and overall productivity of many bird species in the western Sierra Nevada can be relatively high in warm, dry years (Saracco et al., 2019).
Recent technological advances have enabled the rapid expansion of studies identifying geographic structure and tracking movements of populations of small migratory birds across the annual cycle (Mckinnon & Love, 2018;Ruegg et al., 2014;Rushing, Ryder, Scarpignato, et al., 2016); however, a major challenge remains to better understand causes and consequences of movements on population demography, dynamics, and conservation (Ruegg et al., 2020;Saracco & Rubenstein, 2020). Meeting this challenge will be critical for understanding full-life cycle drivers of population change (Hostetler et al., 2015). We suggest that linking individual movement data to demographic monitoring data can be a valuable tool for enhancing our ability to understand drivers of broad-scale dynamics and trends of migratory species. Existing networks of monitoring sites, such as those coordinated by the MAPS program, provide an under-utilized but potentially efficient framework for deploying and recovering tags or obtaining biological samples and providing demographic data for such analyses. The value of these direct measures of movement and demography can be further enhanced by taking advantage of other large-scale observational data sets, such as eBird (Sullivan et al., 2009), which can help to refine understanding of the timing and extent of seasonal movements (Fournier et al., 2017;Vincent et al., 2022).
In addition to the general approach of linking movement and demographic data, our example highlights the potential for leveraging detailed age-specific capture data in a novel integrated population model framework to inform a regional index of average population size, as well as recruitment rates and relative contributions of new yearling recruits vs. older adults and their respective demographic rates in driving population dynamics. Our model builds on previous integrated population models for MAPS data, in which recruitment components of the model were largely latent parameters and for which spatially incongruent data informed the population state (Ahrestani et al., 2017;Saracco & Rubenstein, 2020). As in those previous efforts, years of large population increases appeared to be driven substantially by large recruitment events. Finally, although we focused here on just the survival and recruitment parameters, data on ratios of young to adult birds could also be leveraged to help decompose productivity and juvenile survival components of the recruitment process (Saracco & Rubenstein, 2020).

CO N FLI C T O F I NTE R E S T
None declared.

O PE N R E S E A RCH BA D G E S
This article has been awarded Open Data, Open Materials Badges.
All materials and data are publicly accessible via the Open Science